方程式：

$$ \frac{Du}{Dt}= -\frac{\partial \phi}{\partial x} + \nu\nabla^2 u \\ \frac{Dv}{Dt}= -\frac{\partial \phi}{\partial y} + \nu\nabla^2 v \\ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 $$

渦度：$\displaystyle \zeta \equiv \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ を用いて、渦度方程式： $$ \frac{D\zeta}{Dt}=\nu\nabla^2\zeta $$ が得られる。

非圧縮より、流線関数： $$ u=-\frac{\partial \psi}{\partial y}, v=\frac{\partial \psi}{\partial x} $$ が定義できる。$\zeta=\nabla^2\psi$なので、 $$ \frac{D\nabla^2\psi}{Dt}=\nu\nabla^4 \psi $$ が得られる($\psi$のみの式)。これをスペクトル法で解く（詳しいことは深く説明できないので石岡さんの本を読むべし）。

In [14]:
using FFTW
using Plots
In [2]:
function main(endstep)
    # parameter
    x1 = 4pi/5
    y1 = 4pi/5
    x2 = 6pi/5
    y2 = 6pi/5
    sigma = pi/10
    dt = 1.0/32.0
    freq=32
    N = 256
    K = 170
    I = 2*N
    x = [2*pi*k/I for k=0:I-1]
    #nu = 2.0e-6
    nu = 1.0e-4
    
    
    # Initial condition of vorticity
    zeta0(x,y) = exp( (cos(x-x1)+cos(y-y1)-2)/sigma^2 ) + exp( (cos(x-x2)+cos(y-y2)-2)/sigma^2 ) # I.C. function of zeta
    
    zeta = zeta0.(x', x) # Initial condition of vorticity
    akl_z = fft(zeta)    # Initial condition of vorticity in wavenumber space
    
    # Arrays of wavenumber 
    k1 = [i for i=0:N]; k2 = [i for i=-N+1:-1]
    k_sub = [k1; k2]

    k = zeros(I, I)
    l = zeros(I, I)
    for i = 1:I
        l[i, :] = k_sub
        k[:, i] = k_sub
    end
    
    k2=k.*k
    l2=l.*l
    kl=k.*l
    
    k2l2=k2+l2
    ik2l2=1.0./k2l2
    
    k2_l2=k2-l2
    
    # Initial condition of streamfunction in wavenumber space
    akl = akl_z .*ik2l2
    akl[1,1]=0.0

    # memory allocation
    r=copy(akl)
    s=copy(akl)
    nensei = similar(r)
    uhat = similar(r)
    vhat = similar(r)
    u    =  similar(r)
    v    = similar(r)
    Ahat = similar(r)
    Bhat = similar(r)
    Fhat = similar(r)
    jikan = similar(r)
    k1=similar(r)
    k2=similar(r)
    k3=similar(r)
    k4=similar(r)
    
    # Array for plotting vorticity
    Zetas = Array{Array{Float64,2}}(undef, endstep÷freq+1)
    num=1
    Zetas[num]=Float32.(zeta)

    function RKG!(f, h, y, r, s, args...) # Runge-Kutta-Gill
        r .= h.*f(y, args...)
        y .+= r.*0.5
        s .= h.*f(y, args...)
        y .+= (-1+sqrt(0.5)).*r .+ (1-sqrt(0.5)).*s
        r .= (0.5-sqrt(0.5)).*r .- s
        s .= h.*f(y, args...)
        y .+= r .+ (1+sqrt(0.5)).*s
        r .= -(sqrt(2.0)/3.0*(1.0+sqrt(0.5))).*r .+(2.0/3.0*(-1.0+sqrt(0.5))).*s
        s .= h.*f(y, args...)
        y .+= r .+ s./6.0
    end
    
    function RK4(f, h, x, args...) # 4段4位のRunge-Kutta
        k1 .= h .* f(x,            args...)
        k2 .= h .* f(x + 0.5 * k1, args...)
        k3 .= h .* f(x + 0.5 * k2, args...)
        k4 .= h .* f(x + k3,       args...)
        x = x .+ (k1 .+ 2.0 .* k2 .+ 2.0 .* k3 .+ k4) ./ 6.0
        return copy(x)
    end

    function jikanhatten(akl, nu, k, l) # advection and viscousity
        nensei .= nu.* k2l2 .* akl

        uhat .= -1.0im .* l .* akl
        vhat .=  1.0im .* k .* akl
#         zxhat = -1.0im .* k .* (k.^2 + l.^2) .* akl
#         zyhat = -1.0im .* l .* (k.^2 + l.^2) .* akl

        u .= ifft(uhat)
        v .= ifft(vhat)
#         zx = ifft(zxhat)
#         zy = ifft(zyhat)

#         F = u.*zx .+ v.*zy
#         Fhat = fft(F)
        
        Ahat.=fft(u.*v)
        Bhat.=fft(v.*v.-u.*u)
        Fhat.=-(k2_l2).*Ahat.-kl.*Bhat

        jikan .= Fhat.*ik2l2 .- nensei
        jikan[1,1]=0
        return jikan
    end
    
    
    # main loop
    for itr =1:endstep
        akl[K+1:I-K, :] .= 0.0 + 0.0im
        akl[:, K+1:I-K] .= 0.0 + 0.0im
#         RKG!(jikanhatten, dt, akl, r,s, nu, k, l)
        akl = RK4(jikanhatten, dt, akl, nu, k, l)
        
        if itr%freq==0
            num+=1
            zeta_plot = real.(ifft((k2l2).*akl))
            Zetas[num] = Float32.(zeta_plot)
            print(itr÷freq, " ")
        end

    end
    
    return Zetas
end
np=100
endstep=32*np
@time Zetas=main(1)
@time Zetas=main(endstep);

  2.169812 seconds (5.38 M allocations: 480.841 MiB, 6.22% gc time, 86.28% compilation time)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 553.345757 seconds (2.13 M allocations: 413.989 GiB, 5.38% gc time, 0.06% compilation time)

553秒(約9分)かかった。多分8分くらいまでは速くできそう。

1.5年前に作ったときは900秒だったので、多少は改善してる感（ぜったいfortranのほうが速いと思う）

In [3]:
N=length(Zetas)
anim = @animate for t=1:N
    h=heatmap(Zetas[t],clim=(0,1.0),c=:jet,title="Vortex")
    plot(h, size=(700,600))
end
filename="uzu.gif"
gif(anim, filename, fps=10)

[ Info: Saved animation to C:\Users\*****\Desktop\Julia\Vortex\uzu.gif
Out[3]:
In [13]:
using HDF5
h5open("zeta_0to100s.h5", "w") do file 
    for i=1:101
        name="zeta"*string(i-1)
        write(file, name, Zetas[i])
    end
end

コードの詳細は私まで聞いてください。

2023/05/07